using Printf
using LinearAlgebra
using SparseArrays
For calculating several lowest eigensolutions of large sparse matrix (and other iterative solvers).
using IterativeSolvers
This function emulates speye function which is also found in MATLAB.
speye(N::Int64) = sparse(Matrix(1.0I, N, N));
For plotting (need LaTeX and pgfplots to be installed previously)
using PGFPlotsX
Let's consider again harmonic potential, but now it is in 2d: $$ V(x,y) = \frac{1}{2} m \omega^2 \left( x^2 + y^2 \right) $$
We will use Hartree atomic units: $\hbar = 1$, $m = 1$ and we take $\omega = 1$.
We will compare our numerical eigenvalues with nalytical eigenvalues:
$$ E_{n_{x},n_{y}} = \hbar^2 \omega^2 \left( n_{x} + n_{y} + 1 \right) $$with $n_x$ and $n_y$ start from $0, 1, 2, 3, \ldots$.
Let's first initialize the grid in $x$-direction:
Nx = 5
x_min = -5.0
x_max = 5.0
Lx = x_max - x_min
hx = Lx/(Nx-1)
x = zeros(Float64,Nx)
for i in 1:Nx
x[i] = x_min + (i-1)*hx
end
x
Similarly, we can initialize the grid in $y$-direction:
Ny = 3
y_min = -5.0
y_max = 5.0
Ly = y_max - y_min
hy = Ly/(Ny-1)
y = zeros(Float64,Ny)
for j in 1:Ny
y[j] = y_min + (j-1)*hy
end
y
Now, we will define additional variables to describe our two dimensional grid.
We are essentially build a mapping between 2 indices to 1 indices
Npoints = Nx*Ny
r = zeros(2,Npoints)
ip = 0
idx_ip2xy = zeros(Int64,2,Npoints)
idx_xy2ip = zeros(Int64,Nx,Ny)
for j in 1:Ny
for i in 1:Nx
ip = ip + 1
r[1,ip] = x[i]
r[2,ip] = y[j]
idx_ip2xy[1,ip] = i
idx_ip2xy[2,ip] = j
idx_xy2ip[i,j] = ip
end
end
The grid points are contained in the variable r:
for ip in 1:Npoints
@printf("Point #%3d: [%10.3f,%10.3f]\n", ip, r[1,ip], r[2,ip])
end
function pot_harmonic(x::Float64, y::Float64; ω=1.0)
return 0.5 * ω^2 *( x^2 + y^2 )
end
For visualization purpose, we will use a larger grid size:
Nx = 50
x_min = -5.0
x_max = 5.0
Lx = x_max - x_min
hx = Lx/(Nx-1)
x = zeros(Float64,Nx)
for i in 1:Nx
x[i] = x_min + (i-1)*hx
end
#
Ny = 50
y_min = -5.0
y_max = 5.0
Ly = y_max - y_min
hy = Ly/(Ny-1)
y = zeros(Float64,Ny)
for j in 1:Ny
y[j] = y_min + (j-1)*hy
end
#
Npoints = Nx*Ny
r = zeros(2,Npoints)
ip = 0
idx_ip2xy = zeros(Int64,2,Npoints)
idx_xy2ip = zeros(Int64,Nx,Ny)
for j in 1:Ny
for i in 1:Nx
ip = ip + 1
r[1,ip] = x[i]
r[2,ip] = y[j]
idx_ip2xy[1,ip] = i
idx_ip2xy[2,ip] = j
idx_xy2ip[i,j] = ip
end
end
#
Vpot = zeros(Float64,Npoints)
for ip in 1:Npoints
Vpot[ip] = pot_harmonic(r[1,ip], r[2,ip])
end
@pgf Axis({ height = "10cm", width = "10cm", view=(20,10), colorbar },
Plot3(
{
surf,
},
Coordinates(x, y, reshape(Vpot, Nx, Ny) )
)
)
Let's consider another potential, i.e. Gaussian potential in 2d, centered at $\left( x_c, y_c \right)$:
$$ V_{c}(x,y) = -A\mathrm{e}^{-\alpha\left[ (x-x_c)^2 + (y-y_c)^2 \right]} $$function pot_gaussian( x::Float64, y::Float64; A=1.0, α=1.0, xc=0.0, yc=0.0)
return -A*exp(-α*((x-xc)^2 + (y-yc)^2))
end
x1 = 1.0
y1 = 1.0
x2 = -1.0
y2 = -1.0
Vpot_gaussian = zeros(Float64,Npoints)
for ip in 1:Npoints
Vpot_gaussian[ip] = pot_gaussian(r[1,ip], r[2,ip], xc=x1, yc=y1, α=0.75) +
pot_gaussian(r[1,ip], r[2,ip], xc=x2, yc=y2, α=0.75)
end
@pgf Axis({ height = "10cm", width = "15cm", view=(20,10), colorbar },
Plot3(
{
surf,
},
Coordinates(x, y, reshape(Vpot_gaussian, Nx, Ny) )
)
)
We will start our discussion with a representation of Laplacian operators in 2d, which will be built from one dimensional second derivative operators in one dimension.
We need representation of Laplacian operator $\nabla^2$ in two dimension. We can use Kronecker sum, which in turns can be represented using Kronecker products:
$$ \begin{aligned} L & = D^{(2)}_{x} \oplus D^{(2)}_{y} \\ & = D^{(2)}_{x} \otimes \mathbf{I}_{y} + \mathbf{I}_{x} \otimes D^{(2)}_{y} \end{aligned} $$where $D^{(2)}_{x}$ and $D^{(2)}_{y}$ are representation of second differential operators in $x$ dan $y$ direction and $\mathbf{I}_{x}$ and $\mathbf{I}_{y}$ are identity operators (matrix) in $x$ and $y$ direction.
Let's define again a function to build second differential matrix using 3 point finite difference approximation:
function build_D2_matrix(N::Int64, h::Float64)
mat = zeros(Float64,N,N)
for i = 1:N
mat[i,i] = -2.0
if i != N
mat[i,i+1] = 1.0
mat[i+1,i] = mat[i,i+1]
end
end
return mat/h^2
end
Let's try first with a small grid size:
Nx = 4
x_min = -5.0
x_max = 5.0
Lx = x_max - x_min
hx = Lx/(Nx-1)
x = zeros(Float64,Nx)
for i in 1:Nx
x[i] = x_min + (i-1)*hx
end
#
Ny = 5
y_min = -5.0
y_max = 5.0
Ly = y_max - y_min
hy = Ly/(Ny-1)
y = zeros(Float64,Ny)
for j in 1:Ny
y[j] = y_min + (j-1)*hy
end
#
Npoints = Nx*Ny
r = zeros(2,Npoints)
ip = 0
idx_ip2xy = zeros(Int64,2,Npoints)
idx_xy2ip = zeros(Int64,Nx,Ny)
for j in 1:Ny
for i in 1:Nx
ip = ip + 1
r[1,ip] = x[i]
r[2,ip] = y[j]
idx_ip2xy[1,ip] = i
idx_ip2xy[2,ip] = j
idx_xy2ip[i,j] = ip
end
end
D2x = build_D2_matrix(Nx, hx);
D2y = build_D2_matrix(Ny, hy);
nabla2 = kron(D2x,speye(Ny)) + kron(speye(Nx),D2y);
Matrix(nabla2)
Now, we are ready to solve 2d Schrodinger equation. Like in 1d case, we will convert the differential equation into a matrix eigenvalue problem.
There is one problem that arises, namely the size of Hamiltonian matrix is considerably larger than that of the correnspoding 1d case. We usually don't use the eigh (which uses LAPACK) method in LinearAlgebra method. Even if we can use eigh (assuming we have enough computational resources to use it), we don't usually need all eigenpairs (eigenvalues and eigenvectors) obtained from it. What we need is usually one limited to Nstates (which depend on how many electrons we have) lowest eigenpairs. There are special methods tailored to do such task such as Lanczos, Davidson, conjugate gradient method, etc. In this notebook we will use the so-called LOBPCG (locally-optimal block preconditioned conjugate gradient) method available in the IterativeSolvers package.
Nx = 30
x_min = -5.0
x_max = 5.0
Lx = x_max - x_min
hx = Lx/(Nx-1)
x = zeros(Float64,Nx)
for i in 1:Nx
x[i] = x_min + (i-1)*hx
end
#
Ny = 30
y_min = -5.0
y_max = 5.0
Ly = y_max - y_min
hy = Ly/(Ny-1)
y = zeros(Float64,Ny)
for j in 1:Ny
y[j] = y_min + (j-1)*hy
end
#
Npoints = Nx*Ny
r = zeros(2,Npoints)
ip = 0
idx_ip2xy = zeros(Int64,2,Npoints)
idx_xy2ip = zeros(Int64,Nx,Ny)
for j in 1:Ny
for i in 1:Nx
ip = ip + 1
r[1,ip] = x[i]
r[2,ip] = y[j]
idx_ip2xy[1,ip] = i
idx_ip2xy[2,ip] = j
idx_xy2ip[i,j] = ip
end
end
Vpot = zeros(Float64,Npoints)
for ip in 1:Npoints
Vpot[ip] = pot_harmonic(r[1,ip], r[2,ip])
end
D2x = build_D2_matrix(Nx, hx);
D2y = build_D2_matrix(Ny, hy);
nabla2 = kron(D2x,speye(Ny)) + kron(speye(Nx),D2y);
Ham = -0.5*nabla2 + spdiagm( 0 => Vpot );
Find several lowest (say 5) eigenvalues and eigenvectors using LOBPCG method available in IterativeSolvers:
res = lobpcg(Ham, false, 5, maxiter=500)
res.λ
FD2dGrid struct¶A function to build the grid in 1d:
function init_FD1d_grid( X::Tuple{Float64,Float64}, N::Int64 )
x_min = X[1]
x_max = X[2]
L = x_max - x_min
h = L/(N-1)
x = zeros(Float64,N)
for i = 1:N
x[i] = x_min + (i-1)*h
end
return x, h
end
A structure for FD2dGrid:
struct FD2dGrid
Npoints::Int64
Nx::Int64
Ny::Int64
hx::Float64
hy::Float64
x::Array{Float64,1}
y::Array{Float64,1}
r::Array{Float64,2}
idx_ip2xy::Array{Int64,2}
idx_xy2ip::Array{Int64,2}
end
A constructor of FD2dGrid:
function FD2dGrid(
x_domain::Tuple{Float64,Float64},
Nx::Int64,
y_domain::Tuple{Float64,Float64},
Ny::Int64
)
x, hx = init_FD1d_grid(x_domain, Nx)
y, hy = init_FD1d_grid(y_domain, Ny)
Npoints = Nx*Ny
r = zeros(2,Npoints)
ip = 0
idx_ip2xy = zeros(Int64,2,Npoints)
idx_xy2ip = zeros(Int64,Nx,Ny)
for j in 1:Ny
for i in 1:Nx
ip = ip + 1
r[1,ip] = x[i]
r[2,ip] = y[j]
idx_ip2xy[1,ip] = i
idx_ip2xy[2,ip] = j
idx_xy2ip[i,j] = ip
end
end
return FD2dGrid(Npoints, Nx, Ny, hx, hy, x, y, r, idx_ip2xy, idx_xy2ip)
end
Initialize grid ini two dimensions:
domain_x = (-5.0,5.0)
domain_y = (-5.0,5.0)
Nx = 30
Ny = 30
fdgrid = FD2dGrid( domain_x, Nx, domain_y, Ny );
function build_nabla2_matrix( fdgrid::FD2dGrid )
Nx = fdgrid.Nx
hx = fdgrid.hx
Ny = fdgrid.Ny
hy = fdgrid.hy
D2x = build_D2_matrix(Nx, hx)
D2y = build_D2_matrix(Ny, hy)
∇2 = kron(D2x, speye(Ny)) + kron(speye(Nx), D2y)
return ∇2
end
∇2 = build_nabla2_matrix(fdgrid);
Return the potential, using linear indexing.
function pot_harmonic(fdgrid::FD2dGrid; ω=1.0)
Npoints = fdgrid.Npoints
Vpot = zeros(Npoints)
for i in 1:Npoints
x = fdgrid.r[1,i]
y = fdgrid.r[2,i]
Vpot[i] = 0.5 * ω^2 *( x^2 + y^2 )
end
return Vpot
end
Vpot = pot_harmonic(fdgrid);
Build Hamiltonian matrix:
Ham = -0.5*∇2 + spdiagm( 0 => Vpot );
Use LOBPCG to solve the eigenvalue problems.
res = lobpcg(Ham, false, 5)
The eigenvectors:
X = res.X
This is the eigenvalues:
res.λ
| $n_{x}$ | $n_{y}$ | $E_{n_{x},n_{y}}$ |
|---|---|---|
| 0 | 0 | 1.0 |
| 1 | 0 | 2.0 |
| 0 | 1 | 2.0 |
| 1 | 1 | 3.0 |
| 2 | 0 | 3.0 |
| 0 | 2 | 3.0 |
| 2 | 1 | 4.0 |
Need to use reshape.
@pgf Axis({ height = "10cm", width = "10cm", view=(20,10), colorbar },
Plot3(
{
surf,
},
Coordinates(fdgrid.x, fdgrid.y, reshape(res.X[:,1], (fdgrid.Nx, fdgrid.Ny)) )
)
)
@pgf Axis({ height = "10cm", width = "10cm", view=(20,10), colorbar },
Plot3(
{
surf,
},
Coordinates(fdgrid.x, fdgrid.y, reshape(res.X[:,2], (fdgrid.Nx, fdgrid.Ny)) )
)
)
@pgf Axis({ height = "10cm", width = "10cm", view=(20,10), colorbar },
Plot3(
{
surf,
},
Coordinates(fdgrid.x, fdgrid.y, reshape(res.X[:,3], (fdgrid.Nx, fdgrid.Ny)) )
)
)
@pgf Axis({ height = "10cm", width = "10cm", view=(20,10), colorbar },
Plot3(
{
surf,
},
Coordinates(fdgrid.x, fdgrid.y, reshape(res.X[:,4], (fdgrid.Nx, fdgrid.Ny)) )
)
)
@pgf Axis({ height = "10cm", width = "10cm", view=(20,10), colorbar },
Plot3(
{
surf,
},
Coordinates(fdgrid.x, fdgrid.y, reshape(res.X[:,5], (fdgrid.Nx, fdgrid.Ny)) )
)
)